Counterfactual measure for medical diagnosis

ABSTRACT

A method for providing a computer-implemented medical diagnosis includes receiving an input from a user comprising at least one symptom of the user. The method also includes providing the at least one symptom as an input to a medical model, the medical model being retrieved from memory. The medical model includes a probabilistic graphical model comprising probability distributions and relationships between symptoms and diseases. The method also includes performing inference on the probabilistic graphical model to obtain a prediction of the probability that the user has that disease. The method also includes outputting an indication that the user has a disease from the Bayesian inference, wherein the inference is performed using a counterfactual measure.

CROSS-REFERENCE TO RELATED APPLICATION

The present application claims the benefit of priority under 35 U.S.C. § 119 from U.S. Provisional Patent Application No. 62/812,226, entitled “COUNTERFACTUAL MEASURE FOR MEDICAL DIAGNOSIS,” filed on Feb. 28, 2019, the disclosure of which is hereby incorporated by reference in its entirety for all purposes.

FIELD

Embodiments described herein relate to methods and systems for diagnosis using counterfactual measures.

BACKGROUND

Historically there has been much work towards designing algorithms that can perform medical diagnosis, with models containing hundreds of diseases and symptoms. Recently, there has been renewed hope that these algorithms could eventually surpass the accuracy of human doctors. To achieve this, it is necessary to consider how human doctors perform diagnosis. Given the symptoms presented by the patient and/or relevant knowledge of their background (such as age), a doctor attempts to determine the disease or diseases that are the most likely underlying cause of the symptoms presented. This approach is in contrast to the standard method used for algorithmic diagnosis, which involves ranking diseases by their posterior probabilities given the evidence presented, which is referred to as posterior ranking.

BRIEF DESCRIPTION OF THE FIGURES

FIG. 1 is a schematic of a system in accordance with an embodiment;

FIG. 2 is a diagram of a probabilistic graphical model of the type that can be used with the system of FIG. 1;

FIG. 3 is a diagram of a PGM for explaining the counterfactual diagnosis method in accordance with an embodiment;

FIGS. 4(a)-4(c) are diagrams for explaining a twin network;

FIGS. 5(a)-5(e) are diagrams for explaining a twin network;

FIG. 6 is a plot showing a comparison of results for counterfactual ranking against posterior ranking;

FIG. 7 is a table showing a comparison of results for counterfactual ranking, posterior ranking and doctor ranking;

FIG. 8 is a flow chart depicting a method in accordance with an embodiment;

FIG. 9 is a flow chart showing a method in accordance with a further embodiment and

FIG. 10 is a schematic of a system in accordance with an embodiment.

DETAILED DESCRIPTION

In an embodiment, a method for providing a computer implemented medical diagnosis includes receiving an input from a user comprising at least one symptom of the user. The method also includes providing the at least one symptom as an input to a medical model, the medical model being retrieved from memory. The medical model includes a probabilistic graphical model comprising probability distributions and relationships between symptoms and diseases. The method also includes performing inference on the probabilistic graphical model to obtain a prediction of the probability that the user has that disease. The method also includes outputting an indication that the user has a disease from the Bayesian inference, wherein the inference is performed using a counterfactual measure.

The disclosed system provides an improvement to computer functionality by allowing computer performance of a function not previously performed by a computer. Specifically, the disclosed system provides for a diagnosis system that can provide counterfactual diagnosis that can provide better accuracy in diagnosis. Also, some of the embodiments allows a computationally more efficient diagnosis.

Although, medical diagnosis is exemplified in the description, the method can be applied to any diagnostic problem.

In one embodiment, the counterfactual measure is the expected disablement. In a further embodiment, the probabilistic graphical model is a Noisy-OR model.

In a further embodiment, the counterfactual measure is the expected disablement and the probabilistic graphical model is a Noisy-OR model. This combination allows a closed form solution that only requires knowledge of the observational statistics, as is the case for posterior ranking.

In a yet further embodiment, approximate inference techniques are used, for example, importance sampling. In a yet further embodiment, performing inference comprises using a discriminative model pre-trained to approximate the probabilistic graphical model, the discriminative model being trained using samples from said probabilistic graphical model. Performing inference may also include deriving estimates, from the discriminative model, that the user has a disease. Performing inference may also include performing approximate inference on the probabilistic graphical model to obtain an indication that the user has that disease using the estimate from the discriminative model.

In a further embodiment, the results of predictions of the diseases from the probabilistic graphical model are ranked using the counterfactual measure.

In a further embodiment, the probabilistic graphical model may be a so-called a twin network. The use of a twin network reduces computing counterfactual statements to performing Bayesian inference on an associated causal model.

In a further embodiment, the probabilistic graphical model is a twin network, and the counterfactual measure is the expected disablement. In a further embodiment, the probabilistic graphical model is a twin network, the counterfactual measure is the expected disablement and the probabilistic graphical model is a Noisy-OR model.

In a further embodiment, a system for providing a computer implemented medical diagnosis is provided, the system comprising a processor and a memory, the processor being adapted to receive an input from a user comprising at least one symptom of the user. The processor is also adapted to provide the at least one symptom as an input to a medical model, the medical model being retrieved from memory, the medical model comprising a probabilistic graphical model being stored in memory comprising probability distributions and relationships between symptoms and diseases. The processor is also adapted to perform inference on the probabilistic graphical model to obtain a prediction of the probability that the user has that disease. The processor is also adapted to output an indication that the user has a disease from the Bayesian inference, wherein the inference is performed using a counterfactual measure.

In a further embodiment, a method of forming a twin model that is a graphical representation for computing a counterfactual measure for medical diagnosis is provided. The method comprises receiving a set of data; creating a first graphical representation from the set of data, the first graphical representation comprising a plurality of first nodes indicating symptoms, risk factors and diseases, the first graphical representation indicating a relationship between the symptoms, the risk factors and the diseases; creating a second graphical representation, wherein the second graphical representation is a copy of the first graphical representation, the second graphical representation comprising a plurality of second nodes indicating the symptoms, the risk factors and the diseases; combining the first graphical representation and the second graphical representation to create said twin model; and setting an intervention node of the plurality of second nodes in the twin model to a value so that inference can be performed on the twin model to obtain a prediction of the probability that a user of the model has a specified disease.

This approach massively amortizes the standard inference cost of calculating counterfactuals by abduction, action and prediction, which would otherwise be intractable for large clinical diagnostic models (CDMs).

In a yet further embodiment, each node in the plurality of first nodes is associated with a latent variable, and the method further comprises: linking the first graphical representation and the second graphical representation by sharing the latent variables of each node of the plurality of first nodes with a corresponding node in the plurality of second nodes.

In a yet further embodiment, the first graphical representation and the second graphical representation are probabilistic graphical models.

In a yet further embodiment, the first graphical representation and the second graphical representation are Noisy-OR models.

In a yet further embodiment, the method further comprises: merging a node from the plurality of first nodes with its corresponding node in the plurality of second nodes so as to simplify the twin model.

In a yet further embodiment, the method further comprises: merging a node from the plurality of first nodes with its corresponding node in the plurality of second nodes so as to simplify the twin model.

In a yet further embodiment, after setting the intervention node of the plurality of second nodes in the twin model to a value, the method further comprises: when a node from the plurality of second nodes has identical latent variables and identical parent nodes as its corresponding node in the plurality of first nodes, merging the node with its corresponding node so as to simplify the twin model.

FIG. 1 is a schematic of a diagnostic system. In one embodiment, a user 1 communicates with the system via a mobile phone 3. However, any device could be used, which is capable of communicating information over a computer network, for example, a laptop, tablet computer, information point, fixed computer etc.

The mobile phone 3 will communicate with interface 5. Interface 5 has 2 primary functions, the first function 7 is to take the words uttered by the user and turn them into a form that can be understood by the inference engine 11. The second function 9 is to take the output of the inference engine 11 and to send this back to the user's mobile phone 3.

In some embodiments, Natural Language Processing (NLP) is used in the interface 5. NLP helps computers interpret, understand, and then use every day human language and language patterns. It breaks both speech and text down into shorter components and interprets these more manageable blocks to understand what each individual component means and how it contributes to the overall meaning, linking the occurrence of medical terms to the Knowledge Graph. Through NLP it is possible to transcribe consultations, summarise clinical records and chat with users in a more natural, human way.

However, simply understanding how users express their symptoms and risk factors is not enough to identify and provide reasons about the underlying set of diseases. For this, the inference engine 11 is used. The inference engine is a powerful set of machine learning systems, capable of reasoning on a space of more than hundreds of billions of combinations of symptoms, diseases and risk factors, per second, to suggest possible underlying conditions. The inference engine can provide reasoning efficiently, at scale, to bring healthcare to millions.

In an embodiment, the Knowledge Graph 13 is a large structured medical knowledge base. It captures human knowledge on modern medicine encoded for machines. This is used to allow the above components to speak to each other. The Knowledge Graph keeps track of the meaning behind medical terminology across different medical systems and different languages.

In an embodiment, the patient data is stored using a so-called user graph 15.

Clinical diagnostic models (CDMs) can be probabilistic graphical models (PGMs) that model hundreds of diseases, risk factors (the causes of diseases), and symptoms, as binary nodes.

FIG. 2 is a depiction of a graphical model of the type used in the system of FIG. 1.

The graphical model provides a natural framework for expressing probabilistic relationships between random variables, to facilitate causal modelling and decision-making. In the model of FIG. 2, when applied to diagnosis, D stands for diagnosis, S for symptom and RF for Risk Factor. Three layers: risk factors, diseases and symptoms. Risk factor causes (with some probability) influence other risk factors and diseases, diseases causes (again, with some probability) other diseases and symptoms. There are prior probabilities and conditional marginals that describe the “strength” (probability) of connections. In embodiments, Noisy-OR and Noisy-MAX modelling assumptions are used.

Most widely used is the Noisy-OR assumption, as it closely fits the belief about how diseases develop, and allows large PGMs to be efficiently learned and stored. For a variable Y with parents X₁, . . . ,X_(N), Noisy-OR assumes that Y=1 if any of its parents activate it. A parent X_(i) activates its child if i) the parent is ‘on’, X_(i)=1, and ii) the activation succeeds, with an independent probability of failure λ_(Xi,Y). Further details of the Noisy-OR model are provided below.

In this simplified specific example, the model is used in the field of diagnosis. In the first layer, there are three nodes S₁, S₂ and S₃, in the second layer there are three nodes D₁, D₂ and D₃ and in the third layer, there are three nodes RF₁, RF₂ and RF₃.

In the graphical model of FIG. 2, each arrow indicates a dependency. For example, D₁ depends on RF₁ and RF₂. D₂ depends on RF₂, RF₃ and D₁. Further relationships are possible. In the graphical model shown, each node is only dependent on a node or nodes from a different layer. However, nodes may be dependent on other nodes within the same layer.

In an embodiment, the graphical model of FIG. 2 is a Bayesian Network or probabilistic graphical model “PGM”. In this Bayesian Network (or PGM), the network represents a set of random variables and their conditional dependencies via a directed acyclic graph. Thus, in the network of FIG. 2, given full (or partial) evidence over symptoms S₁, S₂ and S₃ and risk factors RF₁, RF₂ and RF₃ the network can be used to represent the probabilities of various diseases D₁, D₂, and D₃.

Causal knowledge is crucial for effective reasoning in many areas of science and medicine. This is due to the fact that causal relations, unlike correlations, allow one to analyse the consequences of interventions and treatments. When performing a primary diagnosis, a clinician aims to explain a patients presented symptoms by determining the most likely underlying diseases causing them. Despite this, all previous approaches to Machine Learning assisted diagnosis, including both model-based Bayesian and Deep Learning approaches, fail to incorporate causal knowledge in the diagnosis procedure. The embodiments described herein detail a new diagnostic method based on counterfactual inference, which captures the causal aspect of diagnosis overlooked by previous approaches. Whilst previous approaches fail to match human doctors, the counterfactual method achieves expert clinical accuracy.

The embodiments described herein are largely concerned with the PGM. Specifically, methods are presented which instead of using posterior ranking, a counterfactual measure is instead used and inference is performed on this.

Structural causal models (SCMs) can represent variables as functions of their direct causes, along with an unobserved ‘noise’ variable that is responsible for their randomness.

To understand the differences, a background on PGMs will be given:

Definition 1. A structural causal model (SCM) specifies:

-   -   1) a set of independent latent variables U={u₁; : : : ; u_(n)}         distributed according to

${{P(U)} = {\prod\limits_{i}{p\left( u_{i} \right)}}},$

-   -   2) a set of observed variables V={v₁; : : : ; v_(n)};     -   3) a directed acyclic graph G, called the causal structure of         the model, whose nodes are the variables U ∪ V; and     -   4) a collection of functions F={f₁; : : : ; , f_(n)} such that         every observed variable is a function of its parents together         with a latent ‘noise’ term;

v _(i) =f(Pa(v _(i));u _(i));

where Pa_(i) denotes the parent nodes of the i^(th) observed variable from G.

The latent noise term appearing in each f_(i) can be suppressed into Pa(v_(i)) by enforcing the convention that every observed node has an independent latent variable as a parent in G without loss of generality. This convention is adopted throughout the following.

As every observed variable has an independent latent variable as a parent, the distribution over latents induces a distribution over observed variables through the functions F.

That is

${p\left( v_{i} \right)} = {\sum\limits_{{u|v_{i}} = {f_{i}{(u)}}}{{P(u)}.}}$

As P(U) factorises over each element of U, the joint distribution over every node in G factors into conditional distributions for each variable given its parents

${P\left( {v_{1},\ldots \mspace{14mu},v_{n},u_{1},\ldots \mspace{11mu},u_{m}} \right)} = {\prod\limits_{i}{{p\left( v_{i} \middle| {P\; a_{i}} \right)}{\prod\limits_{k}{{p\left( u_{k} \right)}.}}}}$

Hence, once all exogenous latent variables are marginalized over, an SCM describes a Bayesian network. Any Bayesian network can be constructed in this manner from some underlying SCM. Conversely, SCMs can be viewed as a refinement of Bayesian networks to include the details of the underlying generative functions—a redundant feature in Bayesian networks when calculating observational probabilities but vital for calculating counterfactuals.

In order to formally define a counterfactual query, a “do-operator,” an interventional primitive is first defined.

Consider a SCM with functions F. The effect of intervention do(X=x) in this model corresponds to creating a new SCM with functions F_(X=x), formed by deleting from F all functions f_(i) corresponding to members of the set X and replacing them with the set of constant functions X=x. That is, the do-operator forces variables to take certain values, regardless of the original causal mechanism. Probabilities involving the do-operator, such as P(Y=y|do(X=x)) corresponds to evaluating ordinary probabilities in the SCM with functions F_(X), in this case P(Y=y).

Counterfactual queries tell us the probability that certain outcomes would have occurred had some precondition been different. Given some observations

=e, the likelihood is calculated that a different outcome

e′ (counter to the facts of

=e) would have been observed given that some hypothetical intervention do(T=t) had taken place. Here do(T=t) is used to denote the intervention to fix variable T=t as defined by Pearl's calculus of interventions [27]. This counterfactual is summarised by the probability

P(ϵ=e′|ϵ=e,do(T=t))

Counterfactuals allow the quantification of the causal affect of hypothetical interventions on latent variables. For example, it could be asked what is the likelihood that a symptom would not be present, given that the symptom is observed and a disease D, in some unknown state, was treated. This is determined by the counterfactual probability P(S=F|S=T, do(D=F)). If D=T is likely to be causing S=T, i.e. the disease is a good causal explanation of the symptom, then this probability will be high. On the other hand, disease hypotheses that are likely given the symptoms but a poor causal explanation will result in low probabilities.

Definition 2 (Counterfactual). Let X and Y be two subsets of variables in V. The counterfactual sentence “Y would be y (in situation U), had X been x,” is the solution Y=y of the set of equations F_(x), succinctly denoted Y_(x) (U)=y.

As with observed variables in Definition 1, the latent distribution P(U) allows one to define the probabilities of counterfactual statements.

$\begin{matrix} {{P\left( {Y_{x} = y} \right)} = {\sum\limits_{{u|{Y_{x}{(u)}}} = y}{{P(u)}.}}} & (1) \end{matrix}$

It is possible to implement an algorithmic procedure for computing arbitrary counterfactual queries. First, the distribution over latents is updated to account for the observed evidence. Second, the do-operator is applied. Third, the new causal model created by the application of the do-operator in the previous step is combined with the updated latent distribution to compute the counterfactual query.

The above steps can be summarised as follows (where

denotes the set of real evidence):

(1) (abduction) The distribution of the exogenous latent variables P(u) is updated to obtain P(u|

).

(2) (action) Apply the do-operation to the variables in set X, replacing the equations X_(i)=f_(i)(Pa(x_(i)),u_(i)) with X_(i)=x_(i) ∀ X_(i) ε X.

(3) (prediction). Use the modified model to compute the probability of Y=y.

The first step, updating the distribution over latents and storing the resulting probabilities, requires a large amount of computational resources and memory—especially as conditioning on evidence induces correlations between initially uncorrelated latents. Moreover, as this step has to be repeated for every new counterfactual query, such computational resources are continually required. Therefore, the first step, in updating the exogenous latents, is intractable for large models with significant tree-width such as diagnostic networks.

In an embodiment, it is possible to avoid this computationally expensive step by reducing the computing of counterfactual statements to performing Bayesian inference on an associated causal model, known as a twin network.

The twin network method employs two interlinked networks, one representing the real world and the other the counterfactual world being queried. These two networks are identical in structure and share the same latent variables—as the latents are not modified by interventions on the observed variables. The observable variables from the original model are duplicated and labeled distinctly, as they may obtain different values in the real and counterfactual models. After the counterfactual world undergoes the desired intervention, computing a counterfactual statement in the original model is reduced to performing Bayesian inference in the twin network, as the variables in the real and counterfactual worlds can be explicitly compared. The construction of twin networks is discussed further below.

In Pearl, J. 2018. Theoretical impediments to machine learning with seven sparks from the causal revolution. arXiv preprint arXiv:1801.04016), Pearl refers to a “hierarchy” of queries that can be asked about probabilistic models: associative/observational, causal, and counterfactual.

Associative queries lie at the first level of the hierarchy, allowing us only to ask purely observational or “what happened” questions.

This is the domain of Bayesian networks and statistics. Causal queries lie at the second level and concern asking interventional, or “what if”, questions, formalised using the do-operator. At the top of the hierarchy are counterfactual, or “why”, questions. These allow us to ask causal questions about hypothetical interventions, whilst taking into account real evidence—often in apparent contradiction to our desired interventions.

Consider the question “I have not taken an Asprin, A=false. If I had, would my headache H have gone away?”.

This query can be written as p(H_(A=true)=false|A=false) where H_(A=a) is ‘the value of H if A had taken value a’.

Note that the evidence state, A=true, differs from the hypothetical intervention state Â=false—a feature that cannot be represented by any purely causal or observational question. Although we posit a different hypothetical state Â=true to the evidence state A=false, the evidence state is still relevant as it provides information about other variables in our model, such as confounding influences between A and H. It is this feature that separates counterfactual queries from causal queries

Diagnosis—Posterior Ranking

Diagnosis can be performed using a CDM by posterior ranking. Given a patients observed risk-factors R and symptoms S, the posteriors of all model diseases are computed, p(Di=T|R,S), and the most likely diseases returned as a differential. Other factors may be used as well as the posterior likelihood in constructing the final differential, such as disease cruelty, but the posteriors are the key ingredient supplied by the CDM.

Bayesian networks represent arguably the best method for performing algorithmic diagnosis—as they are interpretable and can model a very large number of diseases simultaneously.

The next section will explain how they are used to perform diagnosis via the ranking of disease posteriors given evidence, and give examples of where this approach fails to capture the true diagnostic procedure.

FIG. 3 depicts a typical three-layer Bayesian network used for diagnosis. The top layer is populated by risk factor nodes, such as smoking or a family history of heart disease, representing variables that have a direct causal influence on the diseases. The diseases are found in the middle layer with the bottom layer corresponding to symptoms. For each disease, its children are the symptoms caused by that disease.

The top layer nodes are risk factors, which represent population variables, such as smoking and fitness, that have causal influences on diseases—which constitute the middle layer of nodes. The bottom layer correspond to symptoms, which are caused by underlying diseases. Thick wires highlight back-door paths between diseases D₁, D₂ and D_(K).

FIG. 3 is just an example of one network. Simpler diagnostic networks can been employed, such as Bayesian networks which have no risk factor layer and stipulate no correlations between diseases. Networks with greater topological complexity can also be considered, such as those exhibiting disease-disease connections—allowing for inter-disease causal effects.

For the Diagnosis that is performed on three-layer networks as follows.

A patient supplies evidence which may include known risk factors, symptoms and possibly some known diseases.

These evidence states are denoted as D_(obs); R_(obs); S_(obs). These observations are then applied to the network, and candidate (latent) diseases are ranked by their posteriors. That is, candidate diseases are ranked

={D₁; D₂; : : : ; D_(N)}, such that:

P(D _(i) |D _(obs) ,R _(obs) ,S _(obs))≥P(D _(i+1) |D _(obs) ,R _(obs) ,S _(obs)),∀i

This ranked list of diseases is then used to construct a diagnosis differential, a ranked list of diseases that should be tested for. This can be done in many ways, such as simply choosing the top N most likely diseases, or grading diseases in terms of cruelty to the patient and returning the top N for each cruelty grading. Regardless of how these diseases are ranked, the posterior ranking is the primitive ingredient of such differentials.

This approach correctly identifies the most likely diseases given the evidence. However, this is not necessarily what a doctor does when diagnosing a patient. A doctor seeks to diagnose the disease which they believe to be the underlying cause of the presented evidence. In certain cases, the most likely disease given the evidence does not correspond to the underlying cause of the evidence. If the task of diagnosis is to determine the disease which should be treated in order to reduce or remove the presented symptoms, then outputting the most likely disease will in some cases systematically return the wrong ordering of candidate diseases to treat.

The three-layer CDMs as shown in FIG. 3 have a richer probabilistic structure than the two layer BN2O networks, which model only diseases and symptoms, and make unphysical assumptions such as disease priors being independent. Whilst this additional structure makes these newer models far more accurate than their predecessors, unobserved risk factors result in confounding between diseases, leading to spurious conclusions being drawn if one relies solely on observational statistics (see below).

Confounding: Latent risk factors allow for back-door paths between diseases (e.g. the thick edges in FIG. 3.). This can lead to diseases that have little or nothing to do with the presented symptoms being ranked highly in the differential due to their large posterior correlation with the symptoms through backdoor paths.

In the example DAG depicted in FIG. 3, consider the case where s₂; s₄; s₅=true and all other symptoms are un-evidenced, that is=false. Disease D₂ is a parent of all of the evidenced symptoms and could constitute the most likely cause of the underlying symptoms (depending on the details of the model). Disease D₁ is a parent of S₂ only—it cannot constitute a causal explanation for symptoms S₄ and S₅, whereas disease D_(K) is a parent of none of the evidenced symptoms. D_(K) may be disregarded as it may not be a cause of the symptoms. However its posterior can be very high due to the back-door paths that it shares with D₁ and D₂.

Furthermore, if diseases like D₁ and D₂ are rare diseases such as cancer, with small priors and relatively small posteriors, D_(K) can have a large enough posterior (or even prior) to make it impossible to discern rare but relevant diseases like D₁; D₂ from common, irrelevant diseases like D_(K). Likewise the posterior weight of D₁ derives from a mixture of causative correlations with S₂ and backdoor correlations from sharing multiple risk-factor parents with D₂. Typically, it would be unlikely that D₁ is a good candidate disease for diagnosis, but it is simple to construct scenarios whereby p(D₁|ϵ)≥p(D₂|ϵ).

Counterfactual Ranking

In an embodiment, the counterfactual measure proposed—will be termed expected disablement—and is designed to rank diseases based on how well they explain the symptoms presented. As there are often multiple underlying diseases, it is too stringent to ask that all reported positive symptoms would disappear if it were not for a single disease (for example, if this disease is not a parent of one of these symptoms, the probability of disablement evaluates to 0). On the other hand, disease explanations that can explain as many symptoms as possible are favoured. To accommodate these two characteristics, the proposed measure takes the form of an expectation value.

Definition 3 (Expected disablement). The expected disablement of a disease D_(k) with respect to evidence ϵ is the number of “on” symptoms would be expected to be switched off if there was an intervention to ‘cure’ D_(k), given ϵ

${_{D}\left( {D_{k},ɛ} \right)} = {\sum\limits_{S \Subset S_{+}}{{{S_{+} - S}}{p\left( {S_{D_{k} = 0}^{*} = \left. 0 \middle| ɛ \right.} \right)}}}$

Where S₊={s ε S_(obs)s.t.s=1}, S⁻={s ε S_(obs)s.t.s=0}and ϵ={S_(obs),D_(obs),R_(obs)}. Here for notational connivance a superscript * is added to denote the counterfactual variable.

The expected disablement has three incredibly desirable properties that, mark it out as an important diagnostic measure.

Theorem 1 (Expected disablement decomposition). The expected disablement

_(D)(D_(k); ϵ) can be written as:

$\sum\limits_{S \Subset S_{+}}{{{S_{+} - S}}{p\left( {{D_{D_{k} = 0}^{*} = {\left. 0 \middle| D_{k} \right. = 1}},ɛ} \right)}{p\left( {D_{k} = \left. 1 \middle| ɛ \right.} \right)}}$

By Theorem 1, the expected disablement is an expectation value over the set of positively evidenced symptoms that are children of a given node, weighted by the ‘causal effect’ of D_(k) on its children, and posterior of the disease given the evidence. It thus has the following three desirable properties.

1. Diseases with low posteriors are ranked low

${\lim\limits_{{P{({D_{k} = {1|ɛ}})}}\rightarrow 0}{_{D}\left( {D_{k},ɛ} \right)}} = 0$

2. Diseases that explain many symptoms are ranked higher than diseases that explain fewer symptoms:

${\lim\limits_{{S_{+}\bigcap{{Ch}{(D_{k})}}}->\varnothing}{_{D}\left( {D_{k},ɛ} \right)}} = {0\mspace{11mu} {and}}$ _(D)(D_(k), ɛ) ≤ S₊⋂Ch(D_(k))

3. Diseases with low ‘causal effect’ on the symptoms are ranked low:

P(S_(D_(k) = 0)^(*) = 0|ɛ, D_(k) = 1, ɛ) = 0_(D)(D_(k), ɛ) = 0

-   -   where S={s ε S₊s.t.s ε Ch(D_(k))}S=fs

The Causal approx of expected disablement:

$\sum\limits_{d_{k}}\; {{p\left( d_{k} \middle| ɛ \right)}{{S - S_{+}}}{p\left( S \middle| {{do}\left( {D_{k} = d_{k}} \right)} \right)}}$

-   -   this can be calculated without needing to know the functions

Next Noisy-OR models, a specific class of SCMs for Bernoulli variables that are widely employed as clinical diagnostic models, are defined. The Noisy-OR assumption states that a variable Y is the Boolean OR of its parents X₁,X₂, . . . ,X_(n), where the inclusion or exclusion of each causal parent in the OR function is decided by an independent probability or ‘noise’ term. The standard approach to defining Noisy-OR is to present the conditional independence constraints generated by the Noisy-OR assumption.

${p\left( {{Y = \left. 0 \middle| X_{1} \right.},\ldots \ ,\ X_{n}} \right)} = {\prod\limits_{i = 1}^{n}{p\left( {y = \left. 0 \middle| {{only}\left( {X_{i} = 1} \right)} \right.} \right)}}$

where p(Y=0|only(X_(i)=1)) is the probability that Y=0 conditioned on all of its (endogenous) parents being ‘off’ (X=0) except for X_(i) alone. We denote p(Y=0|only (X_(i)=1))=λ_(Xi,Y) by convention.

The utility of this assumption is that it reduces the number of parameters needed to specify a Noisy-OR network to O(N) where N is the number of directed edges in the network. All that is needed to specify a Noisy-OR network are the single variable marginals p(X_(i)=1) and, for each directed edge Xi→Yj, a single λ_(Xi,Yj). For this reason, Noisy-OR has been a standard assumption in Bayesian diagnostic networks, which are typically large and densely connected and so could not be efficiently learned and stored without additional assumptions on the conditional probabilities. The Noisy-OR assumption for SCMs is now defined.

Definition 4 (Noisy-OR SCM). A Noisy-OR network is an SCM of Bernoulli variables, where for any variable Y with parents Pa(Y)={X₁, : : : ; X_(N)} the following conditions hold.

1. Y is the logical OR of its parents, where for each parent X there is a Bernoulli variable U_(i) whose state determines whether the parent is included in the OR function, or if it is ignored

$\begin{matrix} {y = {\underset{i = 1}{\overset{N}{}}\left( {x_{i}{\overset{\_}{u}}_{i}} \right)}} & \; \end{matrix}$

i.e. Y=1 if any parent is on, x_(i)=1, and is not ignored, u_(i)=0 (u _(l)=1 where ‘bar’ denotes the negation of u_(i)).

2. The probability of ignoring the state of a given parent variable is independent of whether any of the other parents have been ignored

${P\left( {u_{1},u_{2},\ldots \mspace{14mu},u_{N}} \right)} = {\prod\limits_{i = 1}^{N}\; {P\left( u_{i} \right)}}$

3. For every node Y their is a parent ‘leak node’ L_(Y) that is singly connected to Y and is always ‘on’, with a probability of ignoring given by λ_(L) _(Y)

The inclusion of assumption 3 is a convention (a leak node is just a parent L for which p(L=1)=1). The leak node represents the probability that Y=1, even if X_(i)=0 ∀ X_(i) ε Pa(Y). For example, the leak nodes allow the situation that a disease spontaneously occurs to be modelled, even if all risk factors that are modelled are absent, or that a symptom occurs but none of the diseases that are modelled have caused it. Let us briefly consider the physical justification for the Noisy-OR assumption in modelling disease/symptom relationships, given Definition 4.

First, consider the assumption (1), that the generative function is a Boolean OR. This is equivalent to assuming that diseases or risk-factors never ‘destructively interfere’. If D₁ is activating symptom S, and so is D₂, then this joint activation never cancels out to yield S=F. If two diseases are both activating a symptom, then the symptom is activated. Assumption 2 is also well justified in the clinical context. It states that a given disease (if present) has a fixed likelihood of activating a symptom, independent of the presence or absence of any other disease.

The above Noisy-OR assumptions allow great simplification of the model complexity, but also are fairly representative of how diseases/symptoms are expected to respond to their parent risk factors/diseases. For example, a symptom will be present if one parent disease, or another parent disease, is present. Furthermore, it is reasonable to expect that the probability that a symptom does not manifest, even if a parent disease is present, is (typically) independent of whether or not another disease is present.

Noisy-OR models are simple to learn from observational statistics. The latent priors, along with the single-variable priors, constitute all of the parameters needed to specify a Noisy-OR model. The following standard notation is used for the latent variable probabilities.

p(u _(X) _(i,) _(Y) _(j)=0)=λ_(X) _(i,) _(Y) _(j)  (4)

Given the conditional probability distributions for single parent-child pairs, determining these latent priors (and hence fitting a Noisy-OR model) is a computationally simple task.

Theorem 2 (Learning Noisy-OR models). Lambdas can be calculated from pairwise conditional probability tables as

$\begin{matrix} {\lambda_{X_{i},Y_{j}} = \frac{p\left( {Y_{j} = {\left. 0 \middle| X_{i} \right. = 1}} \right)}{p\left( {Y_{j} = {\left. 0 \middle| X_{i} \right. = 0}} \right)}} & (5) \end{matrix}$

As Noisy-OR models can be learned from single and two-variable statistics alone, they are particularly appealing for learning diagnostic models from epidemiological data, which usually comes in the form of double-blind studies involving small sets of variables.

As Noisy-OR models are typically presented as Bayesian networks, the above definition of Noisy-OR is non-standard. Below it is shown that the SCM definition yields the Bayesian network definition, (3).

Theorem 3 (Noisy-OR CPT). The conditional probability distribution of a child Y given its parents {X₁, . . . , X_(n)} and obeying Definition 14 is given by

${p\left( {{Y = {\left. 0 \middle| X_{1} \right. = x_{1}}},\ldots \mspace{14mu},{X_{n} = x_{n}}} \right)} = {\prod\limits_{i = 1}^{n}\; \lambda_{X_{i},Y}^{x_{i}}}$

Proof For Y=0, the negation of y, denoted y, is given by

$\overset{\_}{y} = {{\left( {\underset{i = 1}{\overset{N}{}}\left( {x_{i}{\overset{\_}{u}}_{i}} \right)} \right)} = {\underset{i = 1}{\overset{N}{}}\left( {{\overset{\_}{x}}_{i}u_{i}} \right)}}$

The CPT is calculated from the structural equations by marginalizing over the latents, i.e. all latent states that yield Y=0 are summed over. Equivalently, it is possible to marginalize over all exogenous latent states multiplied by the above Boolean function, which is 1 if the condition Y=0 is met, and 0 otherwise.

$\begin{matrix} {{p\left( {{Y = {\left. 0 \middle| {X\; 1} \right. = x_{1}}},\ldots \mspace{14mu},{X_{n} = x_{n}}} \right)} = {\sum\limits_{u_{Y}}\; {\underset{i = 1}{\overset{N}{}}{\left( {{\overset{\_}{x}}_{i}u_{i}} \right){p\left( u_{Y} \right)}}}}} \\ {= {\sum\limits_{u_{X_{i},Y}}{\prod\limits_{X_{i}}\; {\left( {{\overset{\_}{x}}_{i}u_{i}} \right){\prod\limits_{U_{X_{i},Y}}\; {p\left( u_{X_{i},Y} \right)}}}}}} \\ {= {\prod\limits_{X_{i}}{\sum\limits_{u_{X_{i},Y}}{{p\left( u_{X_{i},Y} \right)}\left( {{\overset{\_}{x}}_{i}{\overset{\_}{u}}_{i}} \right)}}}} \\ {= {\prod\limits_{X_{i}}\left\lbrack {{p\left( {u_{X_{i},Y} = 1} \right)} + {{p\left( {u_{X_{i},Y} = 0} \right)}{\overset{\_}{x}}_{i}}} \right\rbrack}} \\ {= {\prod\limits_{X_{i}}\left\lbrack {\lambda_{X_{i},Y} + {\left( {1 - \lambda_{X_{i},Y}} \right){\overset{\_}{x}}_{i}}} \right\rbrack}} \\ {= {\prod\limits_{X_{i}}\lambda_{X_{i},Y}^{x_{i}}}} \end{matrix}$

This is identical to the Noisy-OR cpt (3).

Next, the expected disablement for Noisy-OR models will be derived. Above, there is a discussion of Pearl's algorithmic procedure for computing the probability of arbitrary counterfactuals.

The first step of this procedure, updating the distribution over latents and storing the resulting probabilities, requires a large amount of computational resources and memory. As this step has to be repeated for every new counterfactual query, the computational overheads can be quite costly. To overcome this, a twin network can be used.

Twin Networks: Twin networks are constructed by taking two copies of the original PGM (e.g. SCM), one to represent real variables and another to represent their counterfactual states. The two copies are linked by sharing exogenous latent variables (in the case of Noisy-OR, these are the noise terms u_(i)) between corresponding nodes. The counterfactual is then computed by applying a do-operation on the counterfactual graph and using belief propagation on the resulting model.

The construction of the twin network is illustrated in FIGS. 4 and 5.

FIG. 4 shows a simple example. FIG. 4(a) shows a structural causal model for two observed variables X, Y and their exogenous latents u_(X), u_(Y). FIG. 4 (b) illustrates the construction of the twin network, with X*,Y* representing the counterfactual states of X, Y. FIG. 4(c) shows the twin network under a hypothetical intervention setting do(X*=F). The counterfactual P(Y*|Y, do(X*=F)) is calculated from (c) using standard belief propagation.

FIG. 5 shows another example of how a twin network is formed. FIG. 5 (a) shows the original PGM/SCM. In this example, the PGM has a directed acyclic graph structure, and comprises a feed-forward, three layer network. The top layer nodes represent risk factors R₁ . . . R_(N). The second layer nodes represent diseases D₁ . . . D_(M). The third layer nodes represent symptoms S₁ . . . S_(L). This model assumes no directed edges between nodes belonging to the same layer.

To construct the twin network, first the SCM as shown in FIG. 5(a) is copied so that there are two networks as shown in FIG. 5(b). In FIG. 5(b) the network on the left will encode the real evidence in a counterfactual query. This network is referred to as the ‘real’ graph.

The network on the right in FIG. 5(b) will take the hypothetical interventions and will yield the hypothetical observations. This network is referred to as the ‘hypothetical’ graph.

Pairs of copied variables are referred to as dual nodes. The nodes on the hypothetical graph denoted X* where X is the corresponding ‘dual’ node on the real graph. For example, in FIG. 5(b) R₁ and R*₁ are dual nodes. Also, D₁ and D*₁ are dual nodes.

The twin network of FIG. 5 (b) is constructed such that each node on the real graph shares its exogenous latent with its dual node, so u*_(Xi)=u_(Xi). These shared exogenous latents are depicted as dashed lines in FIG. 5 (b) to (e).

A disease node is selected in the hypothetical graph to perform the intervention on. In the example shown in FIG. 5 (b), the disease node D*₂ is selected in the hypothetical graph to perform the intervention on.

Blue circles indicate observations, whereas red circles indicate interventions.

Applying the do-operation do(D*=0) to the hypothetical disease node D* cuts any directed edges going into D* and fixes D*=0, as shown in FIG. 5 (d) below.

Once the hypothetical intervention has been applied, it is possible to greatly simplify the twin network graph structure via node merging. In SCM' s a node is fully determined given an instantation of all of its parents and its exogenous latent. Hence, if two nodes have identical exogenous latents and parents, they are copies of each other and can be merged into a single node.

By convention, when dual nodes are merged there is a mapping X*

X (dropping the asterisk). Dual nodes which share no ancestors that have been intervened upon can therefore be merged.

As interventions are not performed on the risk factor nodes, all (R_(i), R*_(i)) are merged (as shown in FIG. 5(d). Note that for the sake of clarity the exogenous latents for risk factors are not depicted in FIG. 5 (c)-(d).

Alternatively, because the risk factor nodes (R_(i), R*_(i)) are not children of the intervention node (D*₂ which is a disease node), the risk factor nodes (R_(i), R*_(i)) can be merged before applying the intervention, as shown in FIG. 5 (c).

Following this, all corresponding real/hypothetical disease nodes that have no intervention ancestors can be merged, as their latents and parents are identical (shown in FIG. 5 (d)).

Finally, any symptoms that are not children of the disease that has been intervened upon can be merged, as again their latent and all of their parent variables are identical. In FIG. 5 (d), only S₁ has any evidence applied to it. As the symptoms are at the bottom of the graph, if they are un-evidenced they can be removed (this is true for all graphical models when doing inference). So only S₁ and its dual S₁* remain.

The resulting twin network is depicted in FIG. 5 (e). In FIG. 5(e), D_(M) is shown as being greyed out. This is because it is not a parent of any of the symptoms. Thus, just like the symptoms that were removed as described above, it is at the bottom of the network and unevidenced. So it is superfluous to the task of inference and hence can be removed. This is the case for all nodes that are not parents of the symptoms.

Counterfactuals are then computed from the posteriors on this twin network. For example in FIG. 5(e), the counterfactual p(S_(1D) ₁ ₌₀ =0|S₁=1) is given by p(S₁*=0|S₁=0, D₂*=0).

This approach massively amortizes the standard inference cost of calculating counterfactuals by abduction, action and prediction, which would otherwise be intractable for large clinical diagnostic models (CDMs).

Noisy-OR Twin Networks

A Noisy-OR SCM (for example that is described above) can be used to derive the conditional probability tables (CPTs) on the twin network. The only CPTs that differ from the original Noisy-OR SCM are those for unmerged dual symptom nodes (i.e., children of the intervention node DK). In one embodiment, the diagnostic twin network is feed-forward and as such the disease layer forms a Markov blanket for the symptoms layer, d-separating dual symptom pairs from each other. As a result, the CPTs for dual symptoms and their parent diseases alone are derived.

Lemma 1: For a given symptom S and its hypothetical dual S*, with parent diseases D and under the hypothetical intervention do(D*k=0), the joint conditional distribution on the twin network is given by

${p\left( {s,\left. s^{*} \middle| D_{\backslash k} \right.,D_{k},{{do}\mspace{14mu} \left( {D_{k}^{*} = 0} \right)}} \right)} = \left\{ \begin{matrix} {p\left( {s = \left. 0 \middle|  \right.} \right)} & {{{if}\mspace{14mu} s} = {s^{*} = 0}} \\ 0 & {{{{if}\mspace{14mu} s} = 0},{s^{*} = 1}} \\ {\left( \frac{1}{\lambda_{D_{k},s} - 1} \right){p\left( {{s = \left. 0 \middle| D_{\backslash k} \right.},{D_{k} = 1}} \right)}{\delta \left( {d_{k} - 1} \right)}} & {{{{if}\mspace{14mu} s} = 1},{s^{*} = 0}} \\ {{p\left( {{s = \left. 1 \middle| _{\backslash k} \right.},{D_{k} = 1}} \right)}/\lambda_{D_{K},s}} & {{{{if}\mspace{14mu} s} = 1},{s^{*} = 1}} \end{matrix} \right.$

where δ(d_(k)−1)=1 if D_(k)=1 else 0, and D=

_(\k) ∪ D_(k), i.e.

_(\k) is the set of diseases that are not intervened on.

The proof for the above Lemma can be found in the Appendix.

When comparing counterfactual and posterior diagnosis, improved diagnostic accuracy should not come at the expense of greater computational cost. To overcome this issue, a closed form expression is now derived for the expected disablement in Noisy-OR diagnostic networks that only requires knowledge of the observational statistics, as is the case for posterior ranking.

Lemma 1 reduces the counterfactual query to an observational query on the original Noisy-OR network, with some simple adjustments to the conditional probabilities. We want to evaluate the expected disablement of the hypothetical intervention do(D*_(k)=0), which is defined as

$_{D_{K},ɛ}:={\sum\limits_{S^{\prime}}\; {{{S_{+} - S_{+}^{\prime}}}{p\left( {\left. S^{\prime} \middle| ɛ \right.,{{do}\left( {D_{k} = 0} \right)}} \right)}}}$

-   -   where the sum enumerates over all possible S′, the hypothetical         symptom states that are consistent with the real evidence         and the intervention do(D_(k)=0). S_(±) denote the         positive/negative symptoms in the evidence set S. The         intervention do(D_(k)=0) can never switch on symptoms, so one         need only enumerate over hypothetical symptoms states where         S′₊⊆S₊ as these are the only counterfactual states with non-zero         weight. From this it also follows that for all s ε S⁻⇒s* ε S′⁻.         Expressing the counterfactual (17) as a query on the twin         network, the positive/negative symptoms on the hypothetical         graph are denoted as S*_(±). The counterfactual posterior         in (17) is represented on the twin network as

p(S ₊ *,S ⁻ *|ϵ,do(D _(k)*=0))=p(S ₊ *,S ⁻ *|S ₊ ,S ⁻,

_(obs) =r,do(D _(K)*=0))

Theorem 4 (Noisy-OR expected disablement). The expected disablement of intervention {circumflex over (D)}_(k)*=0 can be written

$_{{\hat{D}}_{k}} = {\frac{1}{p\left( {S_{+},S_{-}} \right)}{\sum\limits_{X \in S_{+}}\; {\left( {- 1} \right)^{{S_{+}} - {X}}\Lambda_{X,k}{\sum\limits_{s \in X}\; {\frac{1 - \lambda_{sk}}{\lambda_{sk}} \times {P\left( {{S_{-} = 0},{X_{{\hat{D}}_{k = 0}} = 0},{D_{k} = \left. 1 \middle| _{obs} \right.}} \right)}}}}}}$ ${{where}\mspace{14mu} X_{{\hat{D}}_{k = 0}}\mspace{14mu} {means}\mspace{14mu} s_{{\hat{D}}_{k} = 0}{\forall\; {s \in {X\mspace{14mu} {and}\mspace{14mu} \Lambda_{X,k}}}}} = {\prod\limits_{s \in X}\; {\lambda_{k,s}.}}$

Alternatively, the theorem can be expressed as follows:

For the Noisy-OR networks described above, the expected disablement of disease D_(k) is given by

$_{D_{K},ɛ} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)}{\sum\limits_{{C \subseteq S},{{{Ch}{(D_{k})}}}}{{C}{p\left( {{C = 0},{\left. {S_{+}\backslash C} \middle| D_{k} \right. = 1},{_{obs} = r}} \right)}{p\left( {D_{k} = {\left. 1 \middle| _{obs} \right. = r}} \right)}\frac{\prod\limits_{S \in C}\; \left( \lambda_{D_{k},{S - 1}} \right)}{\prod\limits_{S \in S_{+}}\; \left( \lambda_{D_{k},S} \right)}}}}$      where  S_(±) ∖ C  denotes  S⁻ = 0, S₊ ∖ C = 1.

The proof for the above equation can be found in the Appendix.

Note that as all common causes in the three-layer Noisy-OR diagnostic network are observed, the probability P(S⁻=0, X_({circumflex over (D)}) _(k) =0, D_(k)=1|R_(obs)) can be identified in terms of purely observational probabilities via the Backdoor criterion (Pearl, J. 2009. Causality (2nd edition). Cambridge University Press).

Experimentally Comparing Counterfactual and Posterior Diagnosis

Existing algorithms fail to incorporate causal knowledge in the diagnostic process. The counterfactual diagnosis measure captures this causal information.

To evaluate the effectiveness of our counterfactual diagnostic measure, a comparison of the diagnostic precision for the described CDM using posterior ranking and expected disablement ranking of diseases is performed.

The following definition of diagnosis is taken: the determination of the diseases that are most likely to be causing the observed symptoms. We devise a minimal criteria for desirable diagnostic measures to reflect this definition; i) M(Dk)→0 smoothly as P(Dk=T|E)→0 (consistency), ii) M(Dk)=0 for diseases that are not ancestors of any positive symptoms and therefore cannot constitute causal explanations for and observed symptoms (causality), and iii) diseases that explain a greater number of positive symptoms should be ranked highly (simplicity). From this criteria counterfactual diagnostic measure is selected as the expected disablement E(Dk),

${\left( {D_{K},ɛ} \right)}:={\sum\limits_{S^{\prime}}\; {{{S_{+} - S_{+}^{\prime}}}{p\left( {\left. S^{\prime} \middle| ɛ \right.,{{do}\left( {D_{k} = F} \right)}} \right)}}}$

Where E is the total evidence state, S+ denotes the set of positive symptom evidence, S′ is a counterfactual symptom state with S′+={S ε S′|S=T}. The expected disablement is the number of symptoms that one would expect to nullify, had one intervened to cure disease D_(k).

The performance of the diagnostic measure are tested against a set of 600 clinical vignettes. These vignettes are produced by a panel of expert medical professionals, and contain a single ground-truth disease, a collection of symptoms and risk factors, and basic demographic information such as age and gender.

The symptom and risk factor evidence are selected by the expert panel to represent the evidence that could typically be presented to a doctor by a patient who has the vignette disease. Approximately half of the vignettes attempt to model rare diseases, whereas the other half model diseases of high incidence. The task of the diagnostic algorithm is to return a differential of N diseases ranked in terms of their likelihood given the evidence on each vignette. The same diagnostic model is used for evaluating each measure. This model is a Noisy-OR model containing 1300 nodes, constructed using expert knowledge and epidemiological data (Razzaki, S.; Baker, A.; Perov, Y.; Middleton, K.; Baxter, J.; Mullarkey, D.; Sangar, D.; Taliercio, M.; Butt, M.; Majeed, A.; et al. 2018. A comparative study of artificial intelligence and human doctors for the purpose of triage and diagnosis. arXiv preprint arXiv:1806.10698.).

As a simple test of diagnostic accuracy, the accuracy of the top-N differential for posterior and counterfactual ranking were compared, with respect to the ground truth disease. The accuracy is defined as the probability that the ground truth disease is within a differential of size N, for N=1; : : : ; 20.

Given a differential

={D₁,D₂ . . . D_(N)} where

(D_(i))≥

(D_(i+1)) where

(.) is the chosen ranking measure, the accuracy of the ranking is simply calculated as the probability of the inclusion of the ground-truth disease

.

FIG. 6 shows the top-N accuracy of a Noisy-OR test model as calculated on 750 vignettes. As half of vignettes cover rare diseases, for which the ground truth disease should typically not be at the top of the ranking, rankings of N≥5 are focussed upon.

The top line represents the accuracy of differentials of size N for the counterfactual ranking, given by

_(D)(D; ϵ), and likewise for the lower line and posterior ranking P(D; ϵ).

In this regime, using the counterfactual ranking approximately halves the error rate of the diagnosis compared to posterior ranking.

The two diagnostic measures are compared to the diagnoses of human doctors. For each case-card a diagnosis is gathered from four expert doctors, who are given each case card with the model disease masked and asked to provide a list of possible model diseases ranked from very likely to unlikely. All diseases that at least one doctor flagged as likely (P(D=T|

≥0.3) are combined, and this is used as the ‘ground truth’ optimal diagnosis, as it represents the outcome of a 4th opinion expert consultation. For each doctor their precision is measured as the average fraction of the ground truth diseases for each case card that they predicted as likely, and likewise for the clinical diagnostic model (CDM) with M(Dk, E)≥0.3. Whilst this approach has an inherent bias towards doctors, as their assessment define the ground truth, it is necessary in order to model the inherent diagnostic uncertainty. Results are presented in the table of FIG. 7.

In general, the single limiting factor of counterfactual diagnosis is the requirement that the full structural causal model be known in order to compute the expected disablement. Indeed, in purely observational graphical diagnostic models, where the generative functions are unknown, the expected disablement may not be identifiable. However, there is a halfway house between posterior and full counterfactual diagnosis, given by interventional queries from the second level of the causal hierarchy. A simplified version of the quantity envisaged is as follows:

|P(S _(D=0)=0)−P(S=0)|P(D=1|S=1)

As all common causes in a three-layer Noisy-OR diagnostic network are observed, the above quantity can be identified in terms of purely observational probabilities via the Backdoor criterion, assuming the graphical structure of the network is known.

The above counterfactual diagnosis method can be used with any of the known approximate inference techniques, for example, Importance Sampling. The above inference technique can also be used with inference techniques described in Douglas et al “A Universal Marginalizer for Amortized Inference in Generative Models” arXiv:1711.00695v1 [cs.LG] 2 Nov. 2017.

FIG. 8 is a flow chart depicting an embodiment of the counterfactual diagnosis method. In step S801, data is received. This data can be epidemiological data, expert knowledge, data from research studies or other data. A PGM, or SCM, is created (S802). The SCM may be a Noisy-OR SCM. A copy of the SCM is made in S803.

In step S804, the twin network is created using the original SCM from S802 and the copy from S803. Once the twin model is created, the hypothetical intervention can be applied (S805). Next, the twin network can be simplified through merging nodes (S806).

The twin network can then be used in the computation of counterfactuals (S807).

In step S808, the resultant PGM is then used for triage, diagnosis etc. For, example, this can be as described with reference to FIG. 9.

In step S901 of FIG. 9, the user inputs the phrase, for example, “I have a rash” into, for example, the mobile telephone of FIG. 1. This phrase is then passed to system and used to determine a counterfactual query (S902). The counterfactual query is then computed (S903) on the PGM through by identifying a node in the PGM that relates to the user's query. The PGM then outputs the probabilities that each modelled disease is causing the rash (S904).

The system may also determine further questions that will allow the possible causes for the user's symptoms to be narrowed (S905).

While it will be appreciated that the above embodiments are applicable to any computing system, an example computing system is illustrated in FIG. 10, which provides means capable of putting an embodiment, as described herein, into effect. As illustrated, the computing system 500 comprises a processor 501 coupled to a mass storage unit 503 and accessing a working memory 505. As illustrated, a counterfactual diagnosis model 513 is represented as software products stored in working memory 505. However, it will be appreciated that elements of the counterfactual diagnosis model 513, may, for convenience, be stored in the mass storage unit 503. Depending on the use, the counterfactual diagnosis model 513 may be used with a chatbot, to provide a response to a user question that requires the survival analysis model.

Usual procedures for the loading of software into memory and the storage of data in the mass storage unit 503 apply. The processor 501 also accesses, via bus 509, an input/output interface 511 that is configured to receive data from and output data to an external system (e.g. an external network or a user input or output device). The input/output interface 511 may be a single component or may be divided into a separate input interface and a separate output interface.

Thus, execution of the counterfactual diagnosis model 513 by the processor 501 will cause embodiments as described herein to be implemented.

The counterfactual diagnosis model 513 can be embedded in original equipment, or can be provided, as a whole or in part, after manufacture. For instance, the counterfactual diagnosis model 513 can be introduced, as a whole, as a computer program product, which may be in the form of a download, or to be introduced via a computer program storage medium, such as an optical disk. Alternatively, modifications to existing causal discovery model software can be made by an update, or plug-in, to provide features of the above described embodiment.

The computing system 500 may be an end-user system that receives inputs from a user (e.g. via a keyboard) and retrieves a response to a query using counterfactual diagnosis model 513 adapted to produce the user query in a suitable form. Alternatively, the system may be a server that receives input over a network and determines a response. Either way, the use of the counterfactual diagnosis model 513 may be used to determine appropriate responses to user queries, as discussed with regard to FIG. 1.

Implementations of the subject matter and the operations described in this specification can be realized in digital electronic circuitry, or in computer software, firmware, or hardware, including the structures disclosed in this specification and their structural equivalents, or in combinations of one or more of them. Implementations of the subject matter described in this specification can be realized using one or more computer programs, i.e., one or more modules of computer program instructions, encoded on computer storage medium for execution by, or to control the operation of, data processing apparatus. Alternatively, or in addition, the program instructions can be encoded on an artificially-generated propagated signal, e.g., a machine-generated electrical, optical, or electromagnetic signal that is generated to encode information for transmission to suitable receiver apparatus for execution by a data processing apparatus. A computer storage medium (e.g., non-transitory carrier medium) can be, or be included in, a computer-readable storage device, a computer-readable storage substrate, a random or serial access memory array or device, or a combination of one or more of them. Moreover, while a computer storage medium is not a propagated signal, a computer storage medium can be a source or destination of computer program instructions encoded in an artificially-generated propagated signal. The computer storage medium can also be, or be included in, one or more separate physical components or media (e.g., multiple CDs, disks, or other storage devices).

In the above embodiments, a purely directed causal discovery algorithm is converted into one that can also detect latent common causes. The modified algorithm is tested extensively on synthetic and real datasets. In the experiments explained above, the modified algorithm maintained the performance of the original algorithm on the directed datasets and allowed algorithms that originally could not detect latent causes to uncover them in both synthetic and real data.

While certain embodiments have been described, these embodiments have been presented by way of example only, and are not intended to limit the scope of the inventions. Indeed the novel methods and systems described herein may be embodied in a variety of other forms; furthermore, various omissions, substitutions and changes in the form of methods and systems described herein may be made without departing from the spirit of the inventions. The accompanying claims and their equivalents are intended to cover such forms of modifications as would fall within the scope and spirit of the inventions.

Proof of Lemma 1

For a given symptom S and its hypothetical dual S*, with parent diseases

and under the hypothetical intervention do(D_(k)*=0), the joint conditional distribution on the twin network is given by

${p\left( {s,{s^{*}_{\backslash k}},D_{k},{{do}\left( {D_{k}^{*} = 0} \right)}} \right)} = \left\{ \begin{matrix} {{{{p\left( {s = {0}} \right)}\mspace{14mu} {if}\mspace{14mu} s} = {s^{*} = 0}}\mspace{355mu}} \\ {{{{0\mspace{14mu} {if}\mspace{14mu} s} = 0},{s^{*} = 1}}\mspace{439mu}} \\ {{{\left( {\frac{1}{\lambda_{_{k},s}} - 1} \right){p\left( {{s = {0_{\backslash k}}},{D_{k} = 1}} \right)}{\delta \left( {d_{k} - 1} \right)}\mspace{14mu} {if}\mspace{14mu} s} = 1},{s^{*} = 0}} \\ {{{{{p\left( {{s = {1_{\backslash k}}},{D_{k} = 1}} \right)}\text{/}\lambda_{D_{k +}S}\mspace{14mu} {if}\mspace{14mu} s} = 1},{s^{*} = 1}}\mspace{146mu}} \end{matrix} \right.$

where δ(d_(k)−1)=1 if D_(k)=1 else 0, and

=

_(\k) ∪ D_(k), i.e.

_(\k) is the set of diseases that are not intervened on.

First note that for this marginal distribution the intervention do(D_(k)*=0) is equivalent to setting the evidence D_(k)*=0 as we specify the full Markov blanket of (s,s*). Let us denote

_(\k) as the set of parents of (s,s*) not including the intervention node D_(k)* or its dual D_(k). We then need to compute the conditional probability

$\begin{matrix} {{p\left( {s,{s^{*}_{\backslash k}},{D_{k}^{*} = 0},D_{k}} \right)} = {\sum\limits_{n_{0}}{{P\left( u_{s} \right)}{p\left( {{s_{\backslash k}},D_{k},u_{s}} \right)}{p\left( {{s^{*}_{\backslash k}},{D_{k}^{*} = 0},u_{s}} \right)}}}} & (9) \end{matrix}$

The simplest, way to calculate these conditional probabilities is by marginalization, i.e. we circulate p(s=0, s*=0|

. D_(k)*=0, D_(k)), p(s=0|

, D_(k)*=0, D_(k)) and p(s*=0|

, D_(k)*=0, D_(k)) and use these to construct all joint starts. For s_(i)=0, the generative functions are given by

$\begin{matrix} {{p\left( {{s = {0{{Pa}(S)}}},u_{s}} \right)} = {{u_{L}\bigwedge}\underset{X_{i} \in {{Pa}{(S)}}}{}\left( {{\overset{\_}{x}}_{i}\bigvee u_{X_{i},S}} \right)}} & (10) \end{matrix}$

First we compute the joint, state.

${{p\left( {{s = {0 - D_{\backslash k}}},D_{k},u_{s}} \right)}{p\left( {{s^{*} = {0_{\backslash k}}},D_{k}^{*},u_{s}} \right)}} = {{{u_{L}\bigwedge u_{L}}\underset{D_{i} \in _{\backslash k}}{}\left( {u_{D_{k},S}{\overset{\_}{d}}_{i}} \right)\underset{D_{j} \in _{\backslash k}}{}{\left( {u_{D_{j},S}\bigvee{\overset{\_}{d}}_{j}} \right)\bigwedge\left\lbrack {u_{D_{k},S}\bigvee{\overset{\_}{d}}_{k}} \right\rbrack\bigwedge\left\lbrack {u_{D_{k},S}\bigvee{\overset{\_}{d}}_{k}^{*}} \right\rbrack}} = {u_{L}\underset{D_{i} \in _{\backslash k}}{}{\left( {u_{D_{i},S}\bigvee{\overset{\_}{d}}_{i}} \right)\left\lbrack {u_{D_{k},S}\bigvee\left( {{\overset{\_}{d}}_{k}^{*}\bigwedge{\overset{\_}{d}}_{k}} \right)} \right\rbrack}}}$

Where we have used the Boolean identities a ∧a=a and a∨(b∧c)=(a∨b)∧(a∨c). Therefore

${p\left( {{s = 0},{s^{*} = {0_{\backslash k}}},D_{k},D_{k}^{*}} \right)} = {{\sum\limits_{u_{s}}{{P\left( u_{s} \right)}{p\left( {{s = {0_{\backslash k}}},D_{k},u_{s}} \right)}{p\left( {{s^{*} = {0_{\backslash k}}},D_{k}^{*},u_{s}} \right)}}} = {\lambda \; {L_{s}\left\lbrack {{\lambda_{D_{k},S}\left( {d_{k}\bigvee d_{k}^{*}} \right)} + {{\overset{\_}{d}}_{k}\bigwedge{\overset{\_}{d}}_{k}^{*}}} \right\rbrack}{\underset{D_{i} \in _{\backslash k}}{\Pi}\left\lbrack {{\lambda_{D_{i},S}d_{i}} + {\overset{\_}{d}}_{i}} \right\rbrack}}}$

Next, we calculate the single-symptom conditionals

$\begin{matrix} {{p\left( {{s = {0 - D_{\backslash k}}},D_{k}} \right)} = {\sum\limits_{u_{s}}{{P\left( u_{s} \right)}{p\left( {{s = {0_{\backslash k}}},D_{k},u_{s}} \right)}}}} \\ {{= {\sum\limits_{u_{L_{S}}}{{p\left( u_{L_{s}} \right)}u_{L_{s}}\underset{D_{i} \in }{\Pi}{\sum\limits_{u_{D_{i},S}}{{p\left( u_{D_{i},S} \right)}{u_{D_{i},S}\bigvee{\overset{\_}{d}}_{i}}}}}}}} \\ {{= {{p\left( {u_{L_{e}} = 1} \right)}\underset{D_{i} \in }{\Pi}{\sum\limits_{u_{D_{i},S}}\left\lbrack {{p\left( {u_{D_{i},S} = 1} \right)} + {{p\left( {u_{D_{i},S} = 0} \right)}{\overset{\_}{d}}_{i}}} \right\rbrack}}}} \\ {{= {\lambda_{L_{s}}{\underset{D_{i} \in }{\Pi}\left\lbrack {{\lambda_{D_{i},S}d_{i}} + d_{i}} \right\rbrack}\mspace{14mu} 11}}} \end{matrix}$

and similar for p(s*=0|

_(\k), D_(k)*). Note that λr+x=λ^(x). We can now express the joint cpd over dual symptom pairs, using the identities p(s=0, s*=1|X)=p(s=0|X)−p(s=0, s*=0|X), p(s=1, s*=0|X)=p(s*=0|X)−p(s=0, s*=0|X) and p(s=1, s*=1|X)=1−p(s=0|X)−p(s*=0|X)+p(s=0, s*=0|X) for arbitrary X.

${p\left( {s,{s^{*}_{\backslash k}},D_{k},D_{k}^{*}} \right)} = \left\{ \begin{matrix} {{{\lambda_{L_{s}}\lambda_{D_{k},S}^{d_{k}\bigvee d_{k}^{*}}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} {if}\mspace{14mu} s} = {s^{*} = 0}}\mspace{225mu}} \\ {{{{{\lambda_{L_{s}}\left\lbrack {\lambda_{D_{k},S}^{d_{k}} - \lambda_{D_{k},S}^{d_{k}\bigvee d_{k}^{*}}} \right\rbrack}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} {if}\mspace{14mu} s} = 0},{s^{*} = 1}}} \\ {{{{{\lambda_{L_{s}}\left\lbrack {\lambda_{D_{k},S}^{d_{k}^{*}} - \lambda_{D_{k},S}^{d_{k}\bigvee d_{k}^{*}}} \right\rbrack}\underset{D_{i} \in {_{\backslash k}\lambda_{D_{i},S}^{d_{i}}}}{\Pi}\mspace{14mu} {if}\mspace{14mu} s} = 1},{s^{*} = 0}}\mspace{101mu}} \\ {{1 - {{\lambda_{L_{s}}\left\lbrack {\lambda_{D_{k},S}^{d_{k}} + \lambda_{D_{k},S}^{d_{k}^{*}} - \lambda_{D_{k},S}^{d_{k}\bigvee d_{k}^{*}}} \right\rbrack}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} {if}\mspace{14mu} s}} = {s^{*} = 1}} \end{matrix} \right.$

As we are always intervening to switch off diseases D_(k)*=0, then d_(k)∨d_(k) and

$\begin{matrix} {{\lambda_{D_{k},S}^{d_{k}} - \lambda_{D_{k},S}^{d_{k}\bigvee d_{k}^{*}}} = 0} & (12) \end{matrix}$

and therefore p(s=0. s*=1|

_(\k), D_(k), D_(k)*=0)=0 as expected (switching off a disease will never switch on a symptom). This simplifies our expression for the conditional distribution to

${p\left( {s,{s^{*}_{\backslash k}},D_{k},{D_{k}^{*} = 0}} \right)} = \left\{ \begin{matrix} {{{\lambda_{L_{s}}\lambda_{D_{k},S}^{d_{k}}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} {if}\mspace{14mu} s} = {s^{*} = 0}}} \\ {{{{0\mspace{14mu} {if}\mspace{14mu} s} = 0},{s^{*} = 1}}\mspace{259mu}} \\ {{{{{\lambda_{L_{s}}\left\lbrack {1 - \lambda_{D_{k},S}^{d_{k}}} \right\rbrack}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} {if}\mspace{14mu} s} = 1},{s^{*} = 0}}\;} \\ {{{{1 - {\lambda_{L_{s}}\underset{D_{i} \in _{\backslash k}}{\Pi}\lambda_{D_{i},S}^{d_{i}}\mspace{14mu} i\mspace{14mu} s}} = 1},{s^{*} = 1}}} \end{matrix} \right.$

This then simplifies to

${p\left( {s,{s^{*}_{\backslash k}},D_{k},{D_{k}^{*} = 0}} \right)} = \left\{ \begin{matrix} {{{{p\left( {{s = {0_{\backslash k}}},D_{k}} \right)}\mspace{14mu} {if}\mspace{14mu} s} = {s^{*} = 0}}\mspace{284mu}} \\ {{{{0\mspace{14mu} {if}\mspace{14mu} s} = 0},{s^{*} = 1}}\mspace{439mu}} \\ {{{{p\left( {{s = {0_{\backslash k}}},{D_{k} = 0}} \right)} - {{p\left( {{s = {0_{\backslash k}}},D_{k}} \right)}\mspace{14mu} {if}\mspace{14mu} s}} = 1},{s^{*} = 0}} \\ {{{{p\left( {{s = {1_{\backslash k}}},{D_{k} = 0}} \right)}\mspace{14mu} {if}\mspace{14mu} s} = {s^{*} = 1}}\mspace{245mu}} \end{matrix} \right.$

We have arrived at expressions for the cpt's over dual symptoms in terms of cpt's on the real graph, and hence our conterfactual query can be computed on the real graph alone. The third term in simplified conditionals, p(s=0|

_(\k), D_(k)=0)−p(s=0|

_(\k), D_(k)), is identical to zero unless the free D_(K) is in t he state D_(k)=1 in the real graph. Using the definition of noisy-or def eq noisy or to give

$\begin{matrix} {{p\left( {{s = {0_{\backslash k}}},{D_{k} = 0}} \right)} = {\frac{1}{\lambda_{D_{k},S}}{p\left( {{s = {0_{\backslash k}}},{D_{k} = 1}} \right)}}} & (15) \end{matrix}$

we recover

$\begin{matrix} {{{p\left( {{s = \left. 0 \middle| _{\backslash k} \right.},{D_{k} = 0}} \right)} - {p\left( {{s = \left. 0 \middle| _{\backslash k} \right.},D_{k}} \right)}} = {\left( {\frac{1}{\lambda_{D_{k},S}} - 1} \right){p\left( {{s = \left. 0 \middle| _{\backslash k} \right.},{D_{k} = 1}} \right)}{\delta \left( {d_{k} - 1} \right)}}} & (16) \end{matrix}$

where d_(k) is the instantiation of D_(k) on the real graph. The term δ(d_(k)−1) can be thought of as fixing the observation D_(k)=1 on the real graph. Finally, form the definition of the noisy-or cpt we have that

p(s=1|

_(\k) ,D _(k)=0)=p(s=1|

_(\k) ,D _(k)=1)/λ_(D) _(x) ,s  (17)

Proof of Theorem 4 (Noisy-OR Expected Disablement)

For the noisy-or networks described in Appendix A.1-A.4. the expected disablement of disease D_(k) is given by

${E\left( {D_{K},\mathcal{E}} \right)} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)}{\sum\limits_{C \subseteq {S_{+}\bigwedge{{Ch}{(D_{k})}}}}{{}{p\left( {{ = 0},{\left. {S_{\pm}\text{\textbackslash}} \middle| D_{k} \right. = 1},{_{obs} = r}} \right)} {p\left( {D_{k} = {\left. 1 \middle| _{obs} \right. = r}} \right)} \frac{\prod\limits_{S \in }\; \left( {\lambda_{D_{k},S} - 1} \right)}{\prod\limits_{S \in S_{+}}\; \lambda_{D_{N},S}}}}}$

where S_(±) \

denotes S⁻=0,S₊ \

=1. From the above discussion, the non-zero contributions to the expected disablement are

$\begin{matrix} {{E\left( {D_{K},\mathcal{E}} \right)} = {\sum\limits_{ \subseteq S_{+}}\; {{}{p\left( {{S_{-}^{*} = 0},{^{*} = 0},{{S_{+}\text{\textbackslash}} = \left. 1 \middle| S_{+} \right.},S_{-},{_{obs} = r},{{do}\left( {D_{k}^{*} = 0} \right)}} \right)}}}} & (21) \end{matrix}$

Applying Bayes rule, and noting the the real evidence states are not children of the intervention node D*_(K), gives

$\begin{matrix} {{E\left( {D_{K},\mathcal{E}} \right)} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)} {\sum\limits_{ \subseteq S_{+}}{{} {p\left( {{S_{-}^{*} = 0},{^{*} = 0},{{S_{+}\text{\textbackslash}} = {1.\; S_{+}}},{\left. S_{-} \middle| _{obs} \right. = r},{{do}\left( {D_{k}^{*} = 0} \right)}} \right)}}}}} & (22) \end{matrix}$

Furthermore, by applying the inclusion-exclusion principle [58], we can expand this sum over the subsets χ ⊆ S₊ \

to give an expression the right hand side of lil eq where all counterfactual variables are in the zero state

${E\left( {D_{K},\mathcal{E}} \right)} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)}{\sum\limits_{ \subseteq S_{+}}{{}{\sum\limits_{ \subseteq {S_{+}\backslash }}{\left( {- 1} \right)^{{{S_{+}\backslash }} - {}}{p\left( {{S_{-}^{*} = 0},{^{*} = 0},{ = 0},S_{+},{\left. S_{-} \middle| _{obs} \right. = r},{{do}\left( {D_{k}^{*} = 0} \right)}} \right)}}}}}}$

Let us now consider the probabilities Q=p(S⁻*=0,

*=0,S\

*=1,S₊,S⁻|R _(ohs)r,do(D_(k)*=0)). We can express these as marginalizations over the disease layer, which d-separate dual symptom pairs from each other

$\begin{matrix} {Q = {\sum\limits_{d,d_{k}}\; {{p\left( {{_{\backslash k} = d},{D_{k} = {\left. d_{k} \middle| _{obs} \right. = r}}} \right)}{\prod\limits_{S \in S_{-}}\; {{p\left( {{S^{*} = 0},{S = {\left. 0 \middle| _{\backslash k} \right. = d}},{D_{k} = d_{k}},{D_{k}^{*} = 0}} \right)} \times {\prod\limits_{S \in }\; {{p\left( {{S^{*} = 0},{S = {\left. 1 \middle| _{\backslash k} \right. = d}},{D_{k} = d_{k}},{D_{k}^{*} = 0}} \right)}{\prod\limits_{S \in {S_{+}\backslash }}\; {p\left( {{S^{*} = 1},{S = {\left. 1 \middle| _{\backslash k} \right. = d}},{D_{k} = d_{k}},{D_{k}^{*} = 0}} \right)}}}}}}}}} & 24 \end{matrix}$

Applying Lemma 1 simplifies this to

$\begin{matrix} {Q = {\sum\limits_{d}\; {{p\left( {{_{\backslash k} = d},{D_{k} = {\left. d_{k} \middle| _{obs} \right. = r}}} \right)}{\prod\limits_{S \in S_{-}}\; {{p\left( {{S = {\left. 0 \middle| _{\backslash k} \right. = d}},{D_{k} = d_{k}}} \right)}{\prod\limits_{S \in }\; {{p\left( {{S = {\left. 0 \middle| _{\backslash k} \right. = d}},{D_{k} = 1}} \right)}{\delta \left( {d_{k} - 1} \right)} \times {\prod\limits_{S \in {S_{+}\backslash }}\; {{p\left( {{S = {\left. 1 \middle| _{\backslash k} \right. = d}},{D_{k} = 1}} \right)}{\prod\limits_{S \in }\; {\left( {\frac{1}{\lambda_{D_{k},S}} - 1} \right){\prod\limits_{S \in {S_{+}\backslash }}\; \frac{1}{\lambda_{D_{k},S}}}}}}}}}}}}}} & 25 \end{matrix}$

Note that the only Q in lil eq 2 that are not multiplied by a factor of zero have

≠0, and so δ(d_(k)−1) is always present. Marginalizing over d_(k) and simplifying gives

$\begin{matrix} {Q = {{p\left( {{S_{-} = 0},{ = 0},{{S_{+}\backslash } = 0},{D_{k} = {\left. 1 \middle| _{obs} \right. = r}}} \right)}\frac{\prod\limits_{S \in }\; \left( {\lambda_{D_{k},S} - 1} \right)}{\prod\limits_{S \in S_{+}}\; \lambda_{D_{k},s}}}} & (26) \end{matrix}$

Substituting this into lil eq 2 gives

$\begin{matrix} {{E\left( {D_{K},\mathcal{E}} \right)} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)} {\sum\limits_{ \subseteq S_{+}}{{} {p\left( {{S_{-} = 0},{ = 0},{{S_{+}\text{\textbackslash}} = 1},{D_{k} = {\left. 1 \middle| _{obs} \right. = r}}} \right)}\frac{\prod\limits_{S \in }\; \left( {\lambda_{D_{k},S} - 1} \right)}{\prod\limits_{S \in S_{+}}\; \lambda_{D_{k},S}}}}}} & (27) \end{matrix}$

Finally, noting that for any S ε

, S ε Ch(D_(k)), λ_(D) _(k,) _(s) −1=0, therefore the expectation can be reduced to one over all subsets of positively evidenced symptoms that are children of D_(k). The final expression for the expected disablement is

${E\left( {D_{K},\mathcal{E}} \right)} = {\frac{1}{p\left( {S_{+},{\left. S_{-} \middle| _{obs} \right. = r}} \right)} {\sum\limits_{ \subseteq {S_{+}\bigwedge{{Ch}{(D_{k})}}}}{{} {p\left( {{ = 0},{\left. {S_{\pm}\text{\textbackslash}} \middle| D_{k} \right. = 1},{_{obs} = r}} \right)} {p\left( {D_{k} = {\left. 1 \middle| _{obs} \right. = r}} \right)}\frac{\prod\limits_{S \in }\; \left( {\lambda_{D_{k},S} - 1} \right)}{\prod\limits_{S \in S_{+}}\; \lambda_{D_{k},S}}}}}$

where S_(±) \

denotes S⁻=0, S_(+ \)

₌₁ and we have applied the chain rule to D_(k). Before proceeding, we can simplify exp dis significantly.

$\begin{matrix} {{E\left( {D_{K},\mathcal{E}} \right)}:={{\sum\limits_{S^{\prime}}\; {{{S_{+} - S_{+}^{\prime}}}{p\left( {\left. S^{\prime} \middle| \mathcal{E} \right.,{{do}\left( {D_{k} = 0} \right)}} \right)}}} = {\sum\limits_{ \subseteq S_{+}}\; {{}{p\left( {{S_{-} = 0},{ = 0},{{S_{+}\text{\textbackslash}} = \left. 1 \middle| \mathcal{E} \right.},{{do}\left( {D_{k} = 0} \right)}} \right)}}}}} & 29 \end{matrix}$

Next, we can expand the joint state using the inclusion/exclusion principle [58], expressing it as margianls over 0 states

$\begin{matrix} {{{E\left( {D_{K},\mathcal{E}} \right)} = {{\sum\limits_{ \subseteq S_{+}}\; {{}{p\left( {{S_{-} = 0},{ = 0},{{S_{+}\text{\textbackslash}} = \left. 1 \middle| \mathcal{E} \right.},{{do}\left( {D_{k} = 0} \right)}} \right)}}} = {\underset{ \subseteq S_{+}}{\quad\sum}\; {}{\sum\limits_{ \subseteq {S_{+}\backslash }}{\left( {- 1} \right)^{{{S_{+}\backslash }} - {}}{p\left( {{S_{-} = 0},{ = 0},{ = \left. 0 \middle| \mathcal{E} \right.},{{do}\left( {D_{k} = 0} \right)}} \right)}}}}}} \mspace{79mu} {{Next},{{{denote}\mspace{14mu} }:={\bigvee{.}}}}} & 30 \end{matrix}$

We with to count the number of occurrences of the marginal p(S⁻=0,

=0|ϵ, do(D_(k)=0)) in partial sum. Clearly,

varies over all subsets of S₊. For fixed

, valid

that generate

varies over all subsets of

, and for given

, χ is fixed to χ=

\

. Therefore, in partial sum, the coefficient of p(S⁻=0,

=0|ϵ, do(D_(k)=0)) is Σ

_(⊆,)

|

|(−1)^(S) ⁺ ^(\)

^(|−|)

^(\)

^(|), As

⊆ S₊ and χ ⊆

, then |S₊\

|−|

\

|=|S₊|−|

|. Finally, using

$\begin{matrix} {{\sum\limits_{k = 0}^{N}\; {k\begin{pmatrix} N \\ k \end{pmatrix}}} = {N\; 2^{N - 1}}} & (31) \end{matrix}$

we recover

$\begin{matrix} {{{E\left( {D_{K},\mathcal{E}} \right)} = {\sum\limits_{ \subseteq S_{+}}{\left( {- 1} \right)^{{S_{+}} - {}}{A}2^{{} - 1}{p\left( {{S_{-} = 0},{ = \left. 0 \middle| \mathcal{E} \right.},{{do}\left( {D_{k} = 0} \right)}} \right)}}}}} & (32) \end{matrix}$ 

1. A method for providing a computer implemented medical diagnosis, the method comprising: receiving an input from a user comprising at least one symptom of the user; providing the at least one symptom as an input to a medical model, the medical model being retrieved from memory, the medical model comprising: a probabilistic graphical model comprising probability distributions and causal relationships between symptoms and diseases; performing counterfactual inference on the probabilistic graphical model to obtain a prediction of the probability that the user has a disease; and outputting a counterfactual measure determined from the counterfactual inference, wherein the counterfactual measure is the expected number of symptoms that would not be present if a disease was treated, wherein the probabilistic graphical model is a twin network, the twin network formed by: creating a first graphical representation from a set of data, the first graphical representation comprising a plurality of first nodes indicating symptoms, risk factors and diseases, the first graphical representation indicating a relationship between the symptoms, the risk factors and the diseases, the first graphical representation indicating real variables; creating a second graphical representation, wherein the second graphical representation is a copy of the first graphical representation, the second graphical representation comprising a plurality of second nodes indicating the symptoms, the risk factors and the diseases, the second graphical representation indicating counterfactual states of the real variables; and combining the first graphical representation and the second graphical representation by the sharing of exogenous causes to create said twin network.
 2. The method of claim 1, wherein a counterfactual measure is an expected disablement.
 3. The method of claim 1, wherein the probabilistic graphical model is a Noisy-OR model.
 4. The method of claim 1, wherein a counterfactual measure is an expected disablement and the probabilistic graphical model is a Noisy-OR model.
 5. The method of claim 1, wherein performing counterfactual inference comprises using a discriminative model pre-trained to approximate the probabilistic graphical model, the discriminative model being trained using samples from said probabilistic graphical model; deriving estimates, from the discriminative model, that the user has a disease; and performing approximate inference on the probabilistic graphical model to obtain an indication that the user has that disease using the estimate from the discriminative model.
 6. The method of claim 1, wherein the results of predictions of the diseases from the probabilistic graphical model are ranked using a counterfactual measure.
 7. (canceled)
 8. (canceled)
 9. The method of claim 1, wherein the probabilistic graphical model is a twin network, a counterfactual measure is an expected disablement and the probabilistic graphical model is a Noisy-OR model.
 10. A system for providing a computer implemented medical diagnosis, the system comprising a processor and a memory, the processor being adapted to: receive an input from a user comprising at least one symptom of the user; provide the at least one symptom as an input to a medical model, the medical model being retrieved from memory, the medical model comprising: a probabilistic graphical model being stored in memory comprising probability distributions and causal relationships between symptoms and diseases, perform counterfactual inference on the probabilistic graphical model to obtain a prediction of the probability that the user has a disease; and output a counterfactual measure determined from the counterfactual inference, wherein, the counterfactual measure is the expected number of symptoms that would not be present if a disease was treated, wherein the probabilistic graphical model is a twin network, the twin network formed by: creating a first graphical representation from a set of data, the first graphical representation comprising a plurality of first nodes indicating symptoms, risk factors and diseases, the first graphical representation indicating a relationship between the symptoms, the risk factors and the diseases, the first graphical representation indicating real variables; creating a second graphical representation, wherein the second graphical representation is a copy of the first graphical representation, the second graphical representation comprising a plurality of second nodes indicating the symptoms, the risk factors and the diseases, the second graphical representation indicating counterfactual states of the real variables; and combining the first graphical representation and the second graphical representation by the sharing of exogenous causes to create said twin network.
 11. A non-transitory carrier medium comprising computer-readable instructions being adapted to cause a computer to perform: receiving an input from a user comprising at least one symptom of the user; providing the at least one symptom as an input to a medical model, the medical model being retrieved from memory, the medical model comprising: a probabilistic graphical model comprising probability distributions and causal relationships between symptoms and diseases; performing counterfactual inference on the probabilistic graphical model to obtain a prediction of the probability that the user has a disease; and outputting a counterfactual measure determined from the counterfactual inference, wherein, the counterfactual measure is the expected number of symptoms that would not be present if a disease was treated, wherein the probabilistic graphical model is a twin network, the twin network formed by: creating a first graphical representation from a set of data, the first graphical representation comprising a plurality of first nodes indicating symptoms, risk factors and diseases, the first graphical representation indicating a relationship between the symptoms, the risk factors and the diseases, the first graphical representation indicating real variables; creating a second graphical representation, wherein the second graphical representation is a copy of the first graphical representation, the second graphical representation comprising a plurality of second nodes indicating the symptoms, the risk factors and the diseases, the second graphical representation indicating counterfactual states of the real variables; and combining the first graphical representation and the second graphical representation by the sharing of exogenous causes to create said twin network.
 12. The non-transitory carrier medium of claim 11, wherein a counterfactual measure is an expected disablement.
 13. The non-transitory carrier medium of claim 11, wherein the probabilistic graphical model is a Noisy-OR model.
 14. A method of forming a twin model that is a graphical representation for computing a counterfactual measure for medical diagnosis, the method comprising: receiving a set of data; creating a first graphical representation from the set of data, the first graphical representation comprising a plurality of first nodes indicating symptoms, risk factors and diseases, the first graphical representation indicating a relationship between the symptoms, the risk factors and the diseases; creating a second graphical representation, wherein the second graphical representation is a copy of the first graphical representation, the second graphical representation comprising a plurality of second nodes indicating the symptoms, the risk factors and the diseases; combining the first graphical representation and the second graphical representation to create said twin model; and setting an intervention node of the plurality of second nodes in the twin model to a value so that inference can be performed on the twin model to obtain a prediction of the probability that a user of the model has a specified disease.
 15. The method of claim 14 wherein each node in the plurality of first nodes is associated with a latent variable, the method further comprising: linking the first graphical representation and the second graphical representation by sharing the latent variables of each node of the plurality of first nodes with a corresponding node in the plurality of second nodes.
 16. The method of claim 14, wherein the first graphical representation and the second graphical representation are probabilistic graphical models.
 17. The method of claim 14, wherein the first graphical representation and the second graphical representation are Noisy-OR models.
 18. The method of claim 14 further comprising: merging a node from the plurality of first nodes with its corresponding node in the plurality of second nodes so as to simplify the twin model.
 19. The method of claim 15 further comprising: merging a node from the plurality of first nodes with its corresponding node in the plurality of second nodes so as to simplify the twin model.
 20. The method of claim 19, wherein after setting the intervention node of the plurality of second nodes in the twin model to a value, the method further comprises: when a node from the plurality of second nodes has identical latent variables and identical parent nodes as its corresponding node in the plurality of first nodes, merging the node with its corresponding node so as to simplify the twin model.
 21. The method of claim 1, wherein the counterfactual inference further comprises: updating a distribution of exogenous latent variables in the twin network based on the at least one symptom; applying a do-operation to a node indicating the disease in the second graphical representation of the updated twin network; and computing a probability that the user would have the disease based on the updated twin network. 